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TRANSIENT CONFORMAL MAPPING METHOD FOR TWO-DIMENSIONAL 
SOLIDIFICATION OF FLOWING LIQUID ONTO A COLD SURFACE 
by Marvin E. Goldstein and Robert Siegel 
Lewis Research Center 

SUMMARY 

A transient conformal mapping method has been developed to determine the config- 
uration of a frozen region formed during two-dimensional solidification. The transient 
shapes of this region are found by mapping the solidified region into a potential plane and 
then determining the time varying conformal transformations between the potential and 
physical planes. The method is applied to the case where solidification is taking place 
on a cold plate of finite width immersed in flowing liquid having a bulk temperature 
above the freezing point. During the transient, the combined energy resulting from con- 
vection to the interface and latent heat of fusion is conducted through the region from the 
moving interface to the plate. The transient results are compared with a simpler quasi- 
steady solution obtained by letting the growth of the frozen region pass instantaneously 
through a succession of steady- state profile shapes. 


INTRODUCTION 

The present report is concerned with developing a method for solving two- 
dimensional transient solidification problems with a convective boundary condition at the 
moving interface. If a flowing warm liquid is placed in contact with a cold surface which 
is at a temperature below the freezing point of the liquid, a frozen region will grow on 
the plate. The flowing liquid supplies energy by convection to the interface formed be- 
tween the solid and liquid phases. This energy along with the latent heat of fusion re- 
leased during the transient growth must be removed by conduction through the solidified 
region to the cooled plate. There is also internal energy removed as the frozen mate- 
rial is subcooled below its freezing point. The temperature at the solid-liquid interface 
is specified throughout the transient growth by the fact that the interface is within a 
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fraction of a degree of being isothermal at the equilibrium freezing temperature even 
when the interface is moving. 

To obtain a solution to the transient solidification problem, it is necessary to solve 
the heat conduction equation within the solidifying layer in order to determine the growth 
of the solidified region. A shape of this region must be found which will result in an 
isothermal interface and at the same time allow the convective, fusion, and subcooling 
energies to be conducted through the frozen layer to the cooled plate. Since the rate of 
freezing is in general nonuniform along the interface, the temperature derivative (i. e., 
the heat flow) in the solid normal to the interface is an unknown function of position and 
time. This is contrasted with the condition at steady state where the only heat flow at 
the interface is by convection from the liquid. The convection is specified in the prob- 
lem, so in the steady- state case the normal derivative at the interface is known. In the 
transient case the interface boundary condition will be a differential equation relating the 
temperature derivative at the interface to the local rate of freezing. An inherent diffi- 
culty in the solution of freezing problems is the nonlinearity introduced by the fact that 
the interface is moving. As a result, it is not possible to utilize superposition to con- 
struct solutions for time varying conditions. 

In the present problem the energy for subcooling solidified material below the freez- 
ing point will be neglected. This is a common assumption in freezing problems since 
the subcooling energy is usually small compared with the heat of fusion that is liberated 
at the interface and then conducted through the frozen region. For freezing in a flowing 
liquid the convective energy supplied must also be conducted from the interface through 
the frozen region, and this can be a large heat flow compared with the energy for sub- 
cooling. Hence, it is usually the case that practically all of the energy flow in the fro- 
zen region arises from that entering the region at the solid-liquid interface. There are 
only some limited cases, for example, where cryogenic coolants are used, for which the 
subcooling can become important. The solutions in references 1 and 2 show that the sub- 
cooling can be neglected within several percent error if the quantity Cp(tj - t w )/X is 
less than about 1. 

With heat capacity neglected, the heat flow in the solidified region is governed by 
the two-dimensional steady-state heat conduction equation (Laplace equation) which must 
be satisfied at each instant within the frozen region. A numerical solution is difficult 
since the shape of the frozen region is unknown and changing with time. To solve La- 
place's equation in two dimensions, conformal mapping can be utilized. The mapping 
method will be developed here for a transient situation, that is, the mapping functions 
will be allowed to vary with time. The method will be demonstrated by analyzing the 
solidification on a cooled plate of finite width immersed in flowing liquid. The transient 
results will be compared with quasi-steady results obtained by having the frozen region 
pass through a series of instantaneous steady- state frozen configurations during the 
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transient growth. The steady-state profiles which form the basis for the quasi-steady 
solution are the same as those in reference 3 where the steady-state mapping method 
was developed for two-dimensional solidification. 

There has been little analytical work done for two-dimensional solidification. One 
example is reference 4 where an analysis is given for freezing inside a square prism 
with the liquid always at the fusion temperature. The transient frozen-layer configura- 
tions were determined approximately. In this instance the liquid was initially at the 
freezing temperature so there was no convective energy transfer to the frozen interface. 
References 5 to 8 consider steady and transient solidified layers in flow channels. Since 
the channels are of symmetric cross section, being either a tube or gap between paral- 
lel plates, there is only one cross-sectional coordinate. In references 5 to 7, the fro- 
zen layers are two-dimensional as they change in thickness along the channel length. 

The analyses, however, neglect axial conduction within the layers so that the solidifica- 
tion portion of the analysis is locally one-dimensional. 


SYMBOLS 


A 

A n 

a 

b 

^initial 

C n 


E 

h 

% 

J e 


K 

K 


n 


L _ t , t.p 

dimensionless half width of plate, - 

k ‘f-V 

value of A for times before start of transient 
time dependent coefficients in mapping 
half width of plate 

time dependent parameter in mapping 
initial value of b at start of transient 
function of /3 and b defined in eq. (27) 
specific heat of solid 

complete elliptic integral of second kind; quantity defined by eq. (38) 

heat transfer coefficient from flowing liquid to frozen interface 

frozen region in Z -plane 

frozen region in W-plane 

complete elliptic integral of first kind 

definite integrals defined in eq. (53) 

thermal conductivity of solidified material 
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Lq frozen region in £ -plane 

n unit outward normal 

p integer, index in eq. (57) 

Q heat flow rate through frozen layer per unit length of plate 

R c dimensionless position vector, Ar_/a 

F position vector to frozen interface 

S dimensionless frozen-layer-liquid interface 

s frozen-layer-liquid interface 

T dimensionless temperature, (t - t w ) / (t ^ - t w ) 
t temperature 

tj freezing temperature 

t^ liquid temperature 

t surface temperature of cold plate 

U intermediate mapping plane 

W analytic function, cp + ixjy 


X, Y 


x ,y 

z 

“n 

^n 

r 

r 


r^J 

V 




0 


dimensionless coordinates, A—, A^ 

a a 

Cartesian coordinates in physical plane 
dimensionless complex physical plane, X + iY 
time dependent coefficients in mapping eq. (23) 
time dependent coefficients defined in eq. (26) 
frozen region in f2-plane 


k tf " t 

length scale parameter, 

ht r t f 


dimensionless gradient operator, — V 

A 

3T + . 9TV 1 
\ 3X 3Y/ 


dimensionless time, 


h 2 (vM 2 

kpX (t f - t w ) 
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9 time 

X latent heat of fusion 

p density of solidified material 

cp real part of W 

rp imaginary part of W 

S2 intermediate mapping plane 

co argument in SI- plane 

Subscript: 

s on frozen interface 

Superscripts: 
ss steady state 

~ (overbar) complex conjugate 


GENERAL TRANSIENT ANALYSIS 

To help fix the general ideas of the analysis, the following physical problem will be 
considered during the development of the analytical method,, Consider a liquid at con- 
stant temperature t^ flowing over an infinitely long flat plate of width 2a as shown in 
figure 1. Suppose that both vertical sides of the plate are insulated and that the plate 
surface is maintained at a uniform temperature t w which is below the freezing temper- 
ature tj of the liquid. Then a frozen region will grow on the plate until the shape and 
size of the region are such that the heat transferred to the frozen interface by the flow- 
ing liquid is exactly balanced by the heat transferred through the frozen region to the 
plate. If the direction of flow of the liquid is parallel to the long dimension of the plate, 
then the heat transfer coefficient h on the surface of the frozen layer can be assumed 
essentially constant. A procedure will be developed here that predicts the transient 
shape and size of this frozen layer as it grows from some initial state to its final equi- 
librium configuration. 

Before proceeding with the analysis in detail, it is helpful to briefly outline some 
aspects of the general method. In reference 3 a conformal mapping method was devised 
to determine the shapes of steady two-dimensional frozen regions. By applying those 
ideas to the present situation, it is found that the physical coordinates of the frozen re- 
gion can be obtained at each instant in time by carrying out the integral 

Z = / £ dW 
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where £ is a temperature derivative function 


1 = _ 9T + t 3T 
Z dX dY 

and W is the complex potential having negative temperature as the real part and lines 
of constant in the direction of heat flow as the imaginary part, that is, 

W = -T + i^ 

The functions £ and W must be related to each other before the integration for deter- 
mining Z can be performed. 

For freezing with convective heat transfer at the interface, the boundary conditions 
provide information about the shape of the image of the frozen region in the £ and W 
complex planes. For example, the moving interface is always an isotherm at the freez- 
ing temperature, and hence it will be a vertical line in the W-plane. The frozen region 
is then drawn in the £- and W-planes, some of the region boundaries, such as the mov- 
ing interface in the £-plane, being unknown functions of time. Then to relate £ and W, 
both the £ and W regions are mapped into a common fixed intermediate region in a 
plane designated as the f2-plane. The mappings will involve unknown functions of time 
since the region in the fl-plane is chosen to be nontime- varying. 

The integral for determining Z is then carried out in terms of the common vari- 
able 0. Thus 


Z(O) = f 5(0) — (O) dO 

J dfi 

The resulting expression for Z contains the unknown time functions resulting from the 
mappings of the time varying frozen region in the 5“ and W-planes into the fixed region 
in the O-plane. These unknown functions are found by substituting this expression for 
Z (evaluated at the moving boundary) into the heat flow boundary condition at the moving 
interface. 

Each of the previous steps will now be carried out. First consider the boundary 
conditions in detail since these will be needed to represent the frozen region in the tem- 
perature derivative and potential planes. 
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Specification of Physical Boundary Conditions 


Since the problem is two-dimensional, let x and y denote coordinates of an arbi- 
trary point in some fixed cross-sectional plane with the origin of the coordinate system 
as shown in figure 1. Let s denote the two-dimensional free surface of the frozen 
layer in this plane and h denote the unit normal to s directed outward from the frozen 

layer. The position vector of an arbitrary point of s is denoted by r. 

s 



Figure 1. - Two-dimensional solidified layer formed on cold plate. 


At the surface s of the frozen layer, the temperature is constant and equal to the 
freezing temperature t^ of the liquid. The local rate at which heat of fusion is being 
liberated at the freezing interface per unit area is equal to 

- 9 *s 

p \ n (1) 

de 

Heat is being supplied to the frozen interface by convective heat transfer from the warm 
liquid at the rate 


h(t z - t f ) (2) 

per unit interface area. Both of the heat fluxes in expressions (1) and (2) must be bal- 
anced by a conduction heat flux away from the interface given by 

kn • Vt (3) 

7 


L 



ii ii bihi iiimia Bimi 


In view of these considerations, the boundary conditions at the freezing interface can be 
written as 


t(r s , 6) = t f 


kn • Vt - h(tj 


- t f ) = pAn 


ffs 

36 


(4a) 

(4b) 


The remaining boundary conditions are along the plane y = 0 at the surface of the 
cooled plate. For 9 > 0 the plate is maintained at constant temperature. Thus 

t(x, 0, 6 ) = t w -a < x < a (5a) 

For the insulated region on either side of the plate there is no heat flow which implies 

— = 0 y = 0 and x < -a, x > a (5b) 

9y 

It is convenient to introduce the following dimensionless quantities: 


. /VVlS 

Vf- V k 


t - 1. 


T = 


w 


t, - 1 

f w 


X = A 


Y = A 

a 


Rc = A — 
S a 


V = — V 


e _ Ah(t ; - ¥ e . h 2 (h - * t ) ! 


pAa 


kpA (t f - t w ) 
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With these definitions the boundary conditions, equations (4) and (5), have the normalized 
form 


T(R S , 0) = 1 


(6a) 


n • VT - h 



(6b) 


T(X, 0, 0) = 0 -A < X < A 


(7a) 


— = 0 Y = 0 and X < -A, X > A (7b) 

9Y 

Specification of Boundary Conditions in Terms of Complex Quantities 

It has been shown (refs. 1 and 2) that the effect of the heat capacity of the frozen 
layer is negligible in most instances since the subcooling energy is small compared with 
the latent heat removed from the moving interface. Hence, for the purposes of this 
analysis the heat capacity will be neglected; therefore, the heat flow in the frozen region 
is governed by the Laplace equation. Thus the temperature is a harmonic function of po- 
sition within the frozen layer. Let cp denote the harmonic function -T, and let x(/ be 
the harmonic conjugate to cp. The complex variable X + iY will be denoted by Z, and 
the complex variable cp + will be denoted by W. The symbol Iq will be used to 
designate the instantaneous dimensionless region occupied by the frozen material in the 
Z- plane. 

Within the frozen layer the potential function W is a function of time 0 and a 
holomorphic function of position Z. Then at each fixed time, 0 = 0 C , the function 

W(Z, 0 = © c ) for all Z6l Q 

is a holomorphic function of the complex variable Z. In view of this the notation 9W/9Z 
will be used to signify the ordinary derivative dW(Z, 0 = 0 c )/dZ. If the convention is 
adopted whereby the real and imaginary parts of a complex number are identified re- 
spectively with the X and Y components of a vector, then 


VT 


9T . 9T 
9X 9Y 


The Cauchy-Riemann equation dty/dX = 9T/9Y implies 
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Hence, 


3W_ 9T | t d± _ 3T t . 3T 

az ax ax ax ay 


VT = - 


aw 

az 


( 8 ) 


Equation (8) will now be used to express the boundary condition, equation (6b), in terms 
of complex quantities that will be needed in the conformal mapping method. 

Since the surface S of the frozen layer is isothermal, the temperature gradient 
must be normal to it. Since the temperature is increasing in the direction of the outward 
drawn normal, the unit normal vector is given by 


n = 


VT 

VT 


Using equation (8) yields 


aw 


aw 

aw 

az 


az 

az 

aw 


aw 

aw 

az 


az 

az 



aw 

2 az_ 


az 

aw 


aw 

aw 

az 

az 

az 

aw 


(9) 


The first term in equation (6b) can be written as 


n • VT 



VT 


Although 


VT = — + i — 

ax 3Y 

in carrying out the dot product the i is an ordinary unit vector and, hence, the number 
-1 does not appear in front of (3T/aY) 2 . * Then 

*For two vectors written as complex numbers, z^ • = (x^ + iy^) • (x 2 + iy 2 ) 

= x x x 2 + y x y 2 = (le (x x + iy x )(x 2 - iy 2 ) = Re z^. 
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( 10 ) 


aw 
az 

The second term in equation (6b) becomes 


VT = 


VT 


= VT = 


VT 


3R q 

n • - - 

az 

aw 

. . iz 


1 

ii — 

a© 

az 

a© 

s 

az 


aw 



aw 


He 9Z az (11) 

a© aw 


evaluated on the interface. 

Equations (10) and (11) are then inserted into the boundary condition, equation (6b), 
to give 


. /) az az 

i + Re — - — 

a© aw 


az 

aw 


( 12 ) 


Mapping of Frozen Region in Z-Plane into Rectangle in W-Plane 

The boundary conditions on the frozen region provide information from which the 
shape of the region in the potential W-plane can be found. The instantaneous frozen re- 
gion Iq is shown in the physical Z-plane in figure 2. It is clear from the boundary 
conditions (eqs. (6a) and (7a)) specifying the temperatures at the frozen interface and 
plate that, since the real part of W is -T, 



11 



-T(X, Y, ©) = Re W(Z, 0) = -1 

Z e FAB 

(13) 

-T(X, Y, 0) = Re W(Z, 0) = 0 

Z eEDC 

(14) 


where the notation FAB represents the set of points along the boundary FAB. The 
boundary condition (eq. (7b)) together with the Cauchy-Riemann condition 


9T _ d}jy 
3Y 3X 


shows that at any instant of time 

f z e FE 


ax 

Z e CB 


Hence, since ^ is the imaginary part of W, 

9m W(Z, 0) = const 

Z e FE 

(15) 

9m W(Z, 0) = const 

Z e CB 

(16) 


* 
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Equations (13) to (16) show that at each instant of time the region I Q in the physical 
plane maps into a rectangular region Jq in the W-plane as shown in figure 3. The 
height of the rectangle Jq varies with time since the height is related to the heat flow 
through the frozen region (this will be shown in connection with eq. (59)) and must be de- 
termined from the solution to the problem. The symmetry of the rectangle with respect 
to the <p-axis follows from the symmetry of the problem with respect to the imaginary 
axis in the physical Z-plane. 

Mapping of Frozen Region in Z-Plane into £-Plane 

Now consider the complex variable £ defined by 

e=— (i7) 

aw 

Then, it follows from the relation for 3W/3Z immediately preceding equation (8) that 

1 _ 3W _ 9T 3T 

? az ax ay 

which together with the temperature derivative boundary conditions can be used to repre- 
sent the frozen region in the i^-plane. Along the insulated portion of the surface, equa- 
tion (7b) shows that 


9m - = 0 


S 


ZeCB 


for< 


Z e EF 


This is shown in figure 4. It follows from the constant temperature boundary condition 
(eq. (7a)) along the plate that 


— = 0 for -A < X < A 

ax 


Hence, 


Re - = 0 Z e EDC 

S 
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Figure A. - Temperature derivative plane, £ = 


(XL 

\dX 


-1 


The symmetry of the problem shows that 


— = 0 for X = 0 

ax 


Hence, 


/?e - = 0 Z e DA 


Since 


1 

c 





we conclude from these conditions that 
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9m 5 = 0 


/?<?£ = 0 


"\ 


Z e CB 
.Z <= EF 
Z e EDC 
Z SDA 


( 18 ) 


This is indicated in figure 4. Equation (9) and the definition (eq. (17)) show that the out- 
ward drawn normal to the surface FAB is in the -? direction. These considerations 
are sufficient to show that the region Iq in the physical plane must map into the region 
Lq in the £-plane shown in figure 4. The shape of the curve BAF and the length of the 
line BCEF depend on the magnitudes of the temperature derivatives and hence the heat 
flow in the solid at the moving interface. These quantities vary with time and are not 
known at this stage of the solution. 


Integral Relating W and £ to Physical Coordinates 

In order to obtain the frozen layer shape, it is convenient to introduce a parametric 
(intermediate) complex variable £2 which will be used to relate W and £. Consider 
the region r in the £2-plane depicted in figure 5. This region is chosen to not change 
with time and is always bounded by a unit semicircle and the real axis. As will be seen 
later, however, it is necessary to allow the length of the line DE to vary with time. 
Suppose that the mapping £2 — W is known which takes at every instant of time the fixed 
region T in the £2-plane into the time variable region J 0 of the W- plane in the man- 
ner indicated by figures 5 and 3. Also suppose that the mapping £2 — £ which takes T 


A 
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into the variable region Lq in the £ -plane .in the manner indicated by figures 5 and 4 is 
known. The locations of the various points in the physical plane can then be found in 
terms of the parametric variable from the formula obtained by integrating equa- 
tion (17), which is 


Z = / £ df2 + function of 0 (19) 

J an 

In particular by letting £2 be along the semicircle in figure 5, the shape of the freezing 
surface is found at each instant of time. Since W is a known function of S2, and Z is 
known as a function of O from equation (19), it is possible to compute the temperature 
at each point of the physical plane. Thus the solution to the problem can be obtained by 
finding the mappings S2 — W and O — £ at each instant of time and then performing the 
integration in equation (19). 


Determination of the Mapping ft - W 

To determine this mapping, it is convenient to introduce an intermediate variable U 
and to map the rectangular region Jq in the W-plane into the upper half of the U-plane 
in the manner indicated in figure 6. An application of the Schwarz- Christoffel trans- 
formation shows that this mapping is defined by 


aw i 

3U K(VTT7)Vu 2 -b 2 i / ^TT 


( 20 ) 


C B A F E 

H 1 ! I 

-1 -b b 1 

Figure 6. - Intermediate U-plane. 
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The parameter b appearing in this mapping is a function of time. 

The mapping which takes the semicircular region T in the 12- plane into the upper 
half of the U-plane in the manner indicated in the figures is defined by 


b(l + 12 2 ) 

V(1 + 122) 2 - (1 - b 2 )(l - 12 2 ) 2 


12 e T 


( 21 ) 


Substituting equation (21) into equation (20) to eliminate U reveals that the mapping 
which takes T into Jq in the manner indicated by figures 5 and 3 is defined by 


aw = au = 2 

30 30 352 K^iC^ya^V-a-b^a -^) 2 


12 e T (22) 


Determination of the Mapping Q -*•£ 

The mapping which takes the region T into Lq is, of course, holomorphic in the 
interior of r since there are no singularities within the region. An examination of fig- 
ures 4 and 5 indicates that this mapping may be expected to be continuous on the bound- 
aries of r since there are no singularities that occur there. Hence, the mapping has 
an analytic continuation to a region which includes the boundaries of T. Since the func- 
tion 12 — £ is real on the real axis of the 12-plane, it is known from the Schwarz re- 
flection principle that the function 12 — £ can be analytically continued to the region that 
is the mirror image of r with respect to the real axis. Hence, £ can be extended ana- 
lytically to a function that is holomorphic in the interior of a unit circle in the 12-plane 
and is continuous on the boundary of the circle. In view of this, the function 12 — £ has 
a Taylor series expansion about the origin of the 12-plane that converges in a circle that 
contains the closure of T. Since £(12) is real for 12 on the real axis in the 12-plane, it 
follows that the coefficients in the series expansion must be real. Since £(12) is pure 
imaginary for 12 on the imaginary axis in the 12-plane, it follows (since even powers of 
i would be real) that only odd powers of 12 can appear in this expansion. Hence, £ can 
be represented by the convergent series 


£(12,0) = -K 





2n+l 


(23) 
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for all in the closure of r with o; n real for n = 0, 1, 2, 3, . . . . Each coefficient 
o> n in this expansion is, in general, a function of time. Thus the mapping S2 — £ will 
be completely determined if the sequence of real functions of time { a> n } can be found. 

To better understand equation (23), introduce polar coordinates |fi| and for the 
f2- plane by letting 


0 = | f2| e iw 


Then equation (23) becomes 

£(«,©) = -k(Vi - b 2 ^ a n |n| 2n+1 e i(2n+1)aJ 

' ' n=0 

= -k( 1 . 1 - b 2 j ^ a n | n| 2n+ * [cos ( 2 n + l)a> + i sin ( 2 n + l)w] 
' ' n =0 


Since on the solid-liquid interface |fi| =1, this becomes on the interface 


C<«, ©) = 



cos 


(2n + l)w + i sin (2n + l)w] £ e FAB 


Thus the interface of unknown shape in figure 4 has been expressed as a series of cosine 
and sine harmonics. 

The functions a n will be found by requiring that £ satisfy the boundary condition 
(eq. (12)) on the free surface FAB, that is, on the unit semicircle in the fi-plane. By 
introducing the definition of £ given by equation (17), the boundary condition (eq. (12)) 
can be written as 


1 + Re 


5(n, e) 8 ?( n > - 0 ) 


ae 


= ?(«,©)! 


fi G FAB 


(24a) 


Using polar coordinates, on the boundary !2 = e iw , results in equation (24a) becoming 


1 + Re £(e lco , 0) - Z ( - = | £( e iw , 0) 

00 J 


0 < Ctf < TT 


(24b) 
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Integration to Obtain Physical Coordinates in Terms ot Q 


Before the « n can be found by use of equation (24b), it is necessary to calculate 
the complex variable Z in equation (24a) in terms of £2 and the coefficients in the ex- 
pansion (eq. (23)) in order to obtain the 9Z/30 term. K the origin of the coordinate 
system is chosen in the physical plane to be at the point D, then equation (19) shows that 


Z («, 0) - A = 



5 — dS2 

an 


where A is the dimensionless length in figure 2. Substituting equations (22) and (23) 
yields 



n r- 


E a n 


n 


2n+l 


|_V(iE n 2 ) 2 - (i - b 2 )(i - n 2 ) 2 n ~° 


dn + A 


z(n, 0 ) = 


oo 

£•■/ 


sr 


n ^ 
y <*y 


ia 


+ A 


n=0 


Q , (1 + y)* - (1 - b*)(l -yY 


Upon carrying out the indicated integration, it is found that (see the appendix for de- 
tails of the integration) 


Z(S2, 0) - A 


_ ln | Wx(fl) + 


(1 + n 2 ) + (i - b 2 )(i - n 2 ) 


E 

n=0 


°n C n 


oo n 1 oo 

E a n E /3 ^ 2r ‘ b E (25a) 

n=l r=0 n=l 


where 
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(25b) 


X(S2) h (1 + S2 2 ) 2 - (1 - b 2 )(l - O 2 ) 2 
and the are defined recursively by 



and the C are defined in terms of the B n by 
n 'r 



It is necessary to impose the restriction that the distance DC remain constant (and 

equal to A) with time. To accomplish this, notice that point D in figure 6 is the point at 

infinity. This can only occur when the denominator of equation (21) is zero. Now equa- 

2 

ting the denominator to zero gives a quadratic equation for O which is solved to show 
that at point D 



Then the point D maps into the point 
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in the £2 plane, where the + sign is chosen so that D will be in the upper half plane. 
This is also equivalent to 


£2 = 



, \l /2 

yTVbM 

t'Tv/ 


as can be shown by multiplying both the numerator and denominator by 



Since D is at the origin of the physical plane, it follows that 


(28) 



Substituting equations (28) and (29) in equation (25a) yields after some algebraic manipu- 
lation 


A ■ b Z “X - £) c n C n 

n=l n=0 

Substituting this into equation (25a) yields 


(30) 


Z(£2, 0) = 



z vs - A 

InVl - b^ I 


In 


bVx(£2) + (1 + £2 2 ) + (1 - b 2 ) (1 - £2 2 ) 
2 Vl - b 2 


+ 


Z “n Z '? n2r 


n=l r=0 


( 31 ) 


To reduce the double series in equation (31) to a more convenient single sum, notice that 
(see ref. 9, p. 57) 


n-1 


n 


3 n+ lfi2r _ ^ ^ r ^n+l+r^r 

n=l r=0 n=0 r=0 n=0 r=0 


Z “n Z ^ n2r ' Z a n+l Z = Z Z <WC‘" 


OO OO 


oo oo 


^2r 


Z Z “n + UrC 1+r ■ Z « 

r=0 n=0 r=0 p=r+l 


Z v? 


Hence upon defining A r by 


\* Z ^ 

p=r+l 


r = 0, 1, 2, 


(32) 


equation (31) becomes 


bA n - A 

Z(fl, ©) = In 


In Vl - b 2 


bVx(n) + (1 + £2 2 ) + (1 - b 2 )(l - fl 2 ) 


2 Vl - b 2 


QQ 

Vx(n) £ 


Sl 2n A 


n 


n=0 


(33) 


which is the expression for Z that will be used to substitute into the boundary condition 
(eq. (24b)). 


Relation for a n Explicitly in Terms of the A n 

The expression for Z given by equation (33) is in terms of the A n while the ex- 
pression for £ given by equation (23) is in terms of the a n . Hence, before substituting 
equation (23) into boundary condition (24b), it is desirable to express the a n in terms 
of the A n . This can be done through the use of equation (32). To this end, equation (27) 
is used to give 


Z “n C n=“0 C 0 + Z “n C n=Y"b 

n=0 n=l D D 


Z (<S - ofk 

n=l 


(1 - b 2 ) 


Z W - tfkl 

n=l 
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In the definitions of /3 (eq. (26)), the jsJ does not appear, that is, /3j = 0. Then 


I 

n=0 


“n C n = 


Z V 3 ? ' Z Vl + “ b2 ) ( Z Vo " Z Vl) " a 0 

n=l n=2 \n=l n=2 


Now use the definition of A r from equation (32) to give 


OO 

Z “n C n - M“o ' (A ° + A l> ' (1 ' b2)(A 0 " A l>] 
n=0 D 


Then by the additional use of equation (32) to obtain Aq, equation (30) becomes 


A = bA 


o - [“0 - < A 0 + A l> - f 1 - b2 )< A o - A l>] 


which is solved for q?q to obtain 


a 0 = 


(bA 0 - A)b 2 

g + A 0 + Aj + (1 - b^)(A 0 - A x ) 

lnVTZ^ 


Now it follows from definitions (26) and (32) that for r > 1 




p=r+2 


p=r+2 


p=r+2 


From equations (32) and (26) the two summations are found to be 


a u 



1 


E V 3 ? = E “p 4 - “r + l <? +1 = A r - “r+l^r 1 = A r ' “r + l 


p=r+2 p=r+l 


b z (r + 1) 


E V?-l = E V?-l " V?-l - a r+lCl 

p=r+2 p=r 


- Vi -°r±* <W ^ ^ ^ + (1 ~ b2> ] 


b 2 r 


b r(r + 1) 


Substituting into equation (35) gives 


~1 + (1 - b 2 )] (2r + 1) Q , 1 _ [l + (1 - b 2 )~| (2r + 1) 


A r + ~T "V" a r 


' r+1 b 2 r + 1 r+1 b 2 (r + 1) b 2 r + 1 ‘ r + 1 b 2 r 


r fl + (1 - b 2 )1 (2r + 1) 


r + 1 b 2 r(r + l)b 2 F+1 r + 1 * 1 


Cancelling terms and solving for a r gives 

a r = b 2 (r + l)A p+1 + [l + (1 - b 2 )] (2r + l)A r + b 2 rA r _ 1 r = 1, 2, . . . (36) 

Upon combining equations (34) and (36), the ct^'s can be written as 

j=o 

where 



E 



In Vl - b 2 


(38) 
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the symbol ( 2 ) denotes the binomial coefficients whose values are 1, 2, and 1 for 
j = 0, 1, and 2, respectively, and S n q is the Kronecker delta. Equation (37) allows us 
to replace the variables a n in equation (23) by the variables and so, in view of 
equation (33), to completely formulate the boundary condition (eq. (24b)) in terms of the 


Terms in Interface Boundary Condition Needed to Determine the A n and b 


The interface boundary condition as given by equation (24a) or (24b) will provide a 
differential equation which determines the A n . Each of the terms to be substituted into 
this condition will now be derived. 


The term 


Re km. 


In order to insert an expression for dZ /dO into the bound- 


\ 90 / 

ary condition given by equations (24), equation (33) is differentiated with respect to time. 
After collecting terms, this gives 


f 


8Z < n . 9 > = _i_J (bA ° : A)b [<i ♦ a 2 ) - (1 - b 2 )(i - O 2 )] 

30 Vip (1 _ b 2 nn vr7 


bb(l - ft 2 ) 2 A n ^ 2n + X(S2) A n J2 2n > 
n=0 n=0 

J 


(bA Q - A)bb 

• # 
A Q b + bA Q 

In 

bVx(fi) + (1 + ft 2 ) + (1 - b 2 )(l - ft 2 ) 

.(1 - b 2 )(lnVl - b 2 

f InV.-bJ 

2 Vl - b 2 


where the dot denotes differentiation with respect to time. By using equation (38), this 
equation can also be written as 
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8Z(0, Q) 1_ 


Eb 


90 \/x(S 2 ) t(l - b 2 ) j=0 


f 1 - (-D^ 1 - fe2 )] fi2j 


+ b *> E A n E ( 2 )(-D j « 2(n+j) + X) E (j)C 1 - t- 1 ^ 1 - b 2 )] « 2(n+ j> ] 

n=0 j=0 n=0 j=0 


+ E In 


bVx(ft) + (1 + ft 2 ) + (1 - b 2 )(l - ft 2 ) 




(39) 


where it was found by differentiating equation (38) that 


bb(bA Q - A) A Q b + bA Q 


E = u + u 1 

(1 - b 2 )flnVl - b 2 J In Vl - b 5 


(40) 


It follows from equation (25b) that on the interface where ft = e ia> for 0 < co < n 


>i, ,\2 


4 , .\2 


X(e iw ) = (l + e 2iCl> ) - (1 - b 2 )(l - e 2ia, J 


. 2 ico r 2 , 2 X .21 - 2icn,., ,2 . 2 \ 

= 4e cos u> + (1 - b ) sm w = 4e (1 - b sin w) 


Hence, 

Vx(e ia> ) = 2e 1W Vl - b 2 sin 2 a) (41) 

which will be used in equation (39). On the interface the In term in equation (39) be- 
comes 
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b^(l - sin^ go) + cos^ a> + 2b Vl - b^ sin^ co cos u) + (1 - b^)^ sin^ go 
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Combining equations (39), (41), and (42) shows that 


0Z(e lfa) , Q) _ 
00 


Eb 


Vl-b 2 sin 2 CO 1(1 - b2 ) 


Y [l - (-l) j (l - b 2 )]e i(2j_1)co 


oo 2 

bt, Y A n ^(^(-dWM- 1 ^ 

n=0 j=0 


+ E A n ElDf 1 ■ (_1)j(1 ■ b 2 )l e i [ 2 ^ n+ ^"^l a) 
n=0 j=0 



r 

b cos 4- - b 2 sin 2 


u> 


+ i< 


u> + tan 


-(1 - b) sin a) 


b Vl - b 2 sin 2 oo 


+ cos w 


(43) 


Equation (23) shows that the expression for £ on the interface is 


?(e iW , 0) = -K ^1 - b 2 ^ Yj a n e i(2n+1) 


CO 


(44) 


By taking the complex conjugate of equation (44), multiplying by equation (43), and then 
taking the real part, an expression is obtained for the term 


Re 


?(e *" ©) m£h Q > 


00 


■] 
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which will be subsequently substituted into boundary condition (24b) 



(Vl - b 2 

Vl - b 2 sin 2 


r 


• oo 

- EJ r£ 

(1 - b^) n=0 


1 

2 [l-(-l) ;i (l-b 2 )jQ! n cos [2(3 -n- l)wl 
j=0 


0000 2 

+ bb X X A n“k 2 COS [ 2 (l + n ' k - l) w ] 

n=0 k=0 j=0 


00 00 2 

+ 22 A n a k 2 (f) T 1 ‘ " b 2 )lcos 1*2 (n + j - k - l)a>l ► 

n=0 k=0 j=0 J 



The term |g(e ia> , Q)| . - An expression for the values of £ on the interface is given 
by equation (44). However, the absolute value of this expression must be taken. This is 
difficult by ordinary means since it would require squaring the real and imaginary parts 
which are infinite series. As a preliminary step to finding the absolute value notice that, 
if 



OO 


- £ b „ z 


then in view of the Cauchy product rule (ref. 9, p. 162) squaring both sides shows that 


n 


£ a n Z " = £ b „ z ” £ b m zm ■ £ zn £ V.-I 


n=0 


n=0 m=0 


n=0 k=0 
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Hence, equating like powers of Z shows that 


n 

a n = Z Wk " = 0, 1, 2 (46) 

k=0 

This set of equations can be solved recursively to determine the b n in terms of the 
a n . Equation (46) will be applied to find a series expansion for |C(e ia) , 0) | to be used in 
the boundary condition (eq. (24b)). 

To this end, notice that equation (44) implies 




a e i(2n+l)w = e ico/2 



V 5 


i2na> 


and 



Hence, if we let 



(47) 


then 



OO 



n=0 


Using equation (46) shows that the are determined in terms of the a n by 

n 

X ^k^n-k = “n 
k=0 


(48) 


(49) 
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When equation (47) is multiplied by equation (48), all imaginary parts cancel out and the 
following is obtained: 


|C(e iw , ©) | = K(Vl-b 2 ) J £ Vm ei2(n ' m)W 

' / ™_n 


n=0 m=0 


= K 


(Vl-b 2 )^ £ V m cos[2(n-m)a>] 

v 7 n=0 m=0 


(50) 


Equation (50) is the expression for the term which is to be substituted into the right side 
of the boundary condition (eq. (24b)). 


Substitution into Boundary Condition and Integration over Cosine Harmonics 

Equations (45) and (50) are to be substituted into equation (24b). This will give a 
relation containing A and b which are functions of time. The equation will also con- 
tain various sine and cosine harmonics of co. In a way analogous to obtaining the coef- 
ficients in an ordinary Fourier series expansion, the equation will be multiplied through 
by cos 2pco and integrated over the interface. By letting p = 0, 1, 2, 3, . . . , a set 
of simultaneous ordinary differential equations will be obtained for the unknowns A r 
and b. 

Before carrying out this procedure, a few preliminary relations are developed. By 
differentiation it is found that, after considerable simplification, 



r 




9 \ 

co + tan"'*' 

-(1 - b 2 ) sin co 


b COS GO 




/ — 

9(0 

L 

_bVl - b 2 sin 2 co + cos co_ 

J 

Vl - b 2 sin 2 co 


and 


9 

9co 


In 


b cos co + V 1 - b 2 sin 2 co 


b sin co 


Vl - b 2 sin 2 co 


It follows by integrating by parts that for p = 0, 1, 2, 


• • 0 
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Equations (45) and (50) are now substituted into the boundary condition (eq. (24b)). 
The result is multiplied through by cos 2pw for p = 0, 1, 2, 3, . . . and integrated 
between co = 0 and ir/2. Using the results in equation (51) leads to the following infinite 
set of first-order ordinary differential equations for the functions A n and b: 


2tt6. 




Eb 


4K 


(VT7) 


E <•»£[> - - b2 >]< K j-n-l + p ^ K j-„-l-p> 


4(1 - b ) n=0 j=0 


bb 


^n ^ a k Xy (j )^”^^ K j+n-k-l+p + K j+n-k-l-p^ 


n=0 k=0 


j=0 


CO 00 


| I \ E Z (?) I 1 - - » 2 >] ^n+j-k-l+p + K n +j -k-l-p) 


n=0 k=0 j=0 


? £ k -p + k -p] ‘ ' I (1 + V °> E V: 

n=0 


n+p 


(52) 


where for n = 0, 1, 2, 


. the K n are defined by 


/ tt/2 

cos 2na> 

Vl - b^ sin^ O) 


du> 


Equations (40) and (38) show that the derivative of E with respect to time is 

bA Q + A Q b E bb 

Ji = H — 

In Vl - b 2 In Vl - b 2 1 ' fe2 

This is substituted into equation (52) to eliminate E, and the result is 


(53) 
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2jr6 p,0 Eb V ^ 

— ^ — + r L, V 


K 


[fr. 


2) (1 - b ) n=0 


E [ l - - ^j-n-l+p + K j-n-l-p) 

U=o 


2b 


In Vl - b" 


{ 1 

[2(n + p) + 1 


K + ± K _ 

n+p 2(n - p) + 1 n_p J 


OQ OO [ ^ 

+ ^ S S a k] X/ (j )("l)^j+n-k-l+p + K j+n-k-l-j 


n=0 k=0 (J=0 


2 vo r 

1 

J 2 L2(k + p) + 1 

In V 1 - a 

OO OO 

r 2 

E\E“i 

cl E 

n=0 k=0 

U=o 

2 vo b2 r 

l 

J 2 L 2 ( k + p) + 1 

In f 1 - b* 


K, + 


K, 


k+p 2(k - p) + 1 k “ p 


K n+j-k-l-p) 


K. „ + K, „ \ \ = -tt( 1 + 6„ n ) y /L/3. 

k+p 2(k - p) + 1 k - p Jj P ’ £6 ' 


n+p 


(54) 


Now if we set 


H 


(0) _ 2b' 

K p 


In 


k i»p , K k-p 




J o L 2 ( k + P) + 1 2(k - p) + lj 

f 1 - b z 


4? P * K k +P + K k- P 


p = 0, 1, 2, . . . and 
k = 0 , ±1 , ±2 , . . . 


> 


(55) 


and 
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min(2, r) 


J £k,p< b >= X f ^ ,1 )<- 1)i [<- 1 > ri -< 1 - b2 >] 

j=o J 


u(l) + 6 H (°) 

ii+j-k-1, p °n, 0 H k, p 


r = 1, 2, 3 and n, k, p = 0, . 1, 2, . . . (56) 


then equation (54) can be put in the form 



OO 

+ 1 V A j( 2 > 

+ b L n n, k i 
D n=0 


n=0 k=0 


t (3) 

a k J n, k, p 


= -’< 1 + 6 p,0 ) Z V 


n=0 


n+p 



P = 0, 1, 2, . . . (57) 


This infinite set of first-order differential equations completely determines the unknown 
functions of time, b and A n , for n = 0, 1, 2, 3, . . . and so in view of equations (33) 
and (22) completely determines the solution to the problem. 


Coordinates of the Solid-Liquid Interface 

If equations (41) and (42) are substituted in equation (33) and the real and imaginary 
parts of the resulting expression are taken, the following parametric equations for the 
shape of the freezing surface are obtained: 
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OO *> 

s 1 - sin^ co ^ cos (2n + l)co 

n=0 


bA Q - A 


In 




b cos 


a* + Vl - b^ sin^ CO 




Y g = 2 T1 - b 2 sin^ ui Y A fi sin (2n + l)w 


n=0 


0 < u> < 

2 

>- (58) 

- — < tan - * < 0 ' 

2 


bAg - A 


In 




co + tan' 


-1 


-(I - b z ) sin a? 


^b Vl - b 2 sin 2 


CO + COS CO/ 




The evaluation of the set of differential equations to obtain the A n and b as func- 
tions of time and the evaluation of the transient interface shapes will be given subse- 
quently. 


Heat Flow Through Frozen Region 

An expression for the heat flow per unit plate length through the frozen layer can be 
obtained from the temperature gradient at the cold plate, 


Q = f k — dx 
•^late 3y 

or in dimensionless form 


Q 

2k(t f - t w ) 



dX 


It has been pointed out in connection with equation (8) that 3T/3Y = 3i^/3X. Hence, 


Q 

2k(t f - t w ) 
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f 

It can be seen from figure 3 that = i/'g - Using the coordinates in the U- 

plane (see fig. 6) gives 

/ -b 

— dU 

A „ 9U 


Substituting equation (20) for 3W/3U gives 



Some special cases will now be considered. These will aid in the interpretation of 
the transient results to be obtained later. 

STEADY-STATE SOLUTION 

The steady-state frozen layer shapes for the geometry considered here have been 
already given in reference 3. They can be obtained from the present analysis by letting 
all the time dependent coefficients A n be zero. Then from equation (58) 



The steady-state profiles are thus a function of A and b. The parameter A contains 
the physical quantities governing the frozen configuration 
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A 


ha t% l " l t 


k tf-*w 


and hence a relation is needed between the mapping parameter b and A. To obtain this 

• • 

relation, note that at steady state, b and A are zero so equation (57) reduces to 


+ <*p, 0> £ V: 


n=0 


n+p 



With all the A r = 0, equation (37) gives 


G '0 = 


Ab 

In Vl - b^ 


o> n = 0 for n > 1 


(61) 


(62) 


Then from equation (49) 


% = “0 /2 - 1/ pLr 4 = 0 for ” 2 1 

1 1„|W 

Equation (61) reduces to the following single relation for p = 0 

_ 2Ab 2 = 0 

In Vl - b^ k(Vi - b 2 j 

so that the relation between b and A at steady state is 


(63) 


A = - 


In Vl - b 2 



(64) 


This is the same as the result obtained in reference 3. This relation has been plotted in 
figure 7. For a given A as determined by the imposed temperatures and heat transfer 
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Mapping parameter, b 
(a) Range of b, 0. 4 to 0. 96. 



J I I L I 1 I I I 

.91 .92 .93 .94 .95 .96 .97 .98 .99 

Mapping parameter, b 

(b) Range of b, 0. 90 to 1.00. 

Figure 7. - Steady-state relation between cooling parameter A and mapping 
parameter b (eq. (64)). 





Mapping parameter, b 
(0 Range of b, 0.990 to 1.000. 

Figure 7. - Concluded. 

coefficient, b can then be found. Equations (60) are evaluated for various values of b, 
and these contours are given in figure 8. Hence, by using figures 7 and 8 the steady- 
state contour of the frozen region can be quickly found for a given A. 

It follows from equation (59) that the heat flow through the layer is only a function of 
b. In fact, for both transient and steady conditions 

§ = KQ» 

2k <‘l-‘w> K / 

This relation between Q/2k(t f - t w ) and b is plotted in figure 9. 

To summarize for the steady-state solution: for given physical conditions compute 
the value of A, and then determine b from figure 7. Figures 8 and 9 then give the con- 
tour of the frozen region and the heat flow through it. 
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Dimensionless coordinate, y la 






Mapping parameter, b 
(a) Range of b, 0.4 to 1.0. 



... .1 l . I _ l. I I . J I i 

.91 .92 .93 .94 .95 .96 .97 .98 .99 

Mapping parameter, b 

(b) Range of b, 0.90 to 1.00. 

Figure 9. - Dimensionless heat flow through frozen layer as a function of mapping parameter b (eq. (65)). 



Figure 9. - Concluded. 


QUASI-STEADY SOLUTION 
Analysis 

If all the A n are zero, then the only o- n that is not zero is a Q . In this case, the 
mapping function (eq. (23)) reduces to 


5 = -K(Vl - b 2 )a 0 ft 


On the frozen interface where = e 


1 CO 


£ = -K 




«Qe iaj 0 — co < 7r 


This is a semicircle which is the shape of the mapping for steady state. Thus, if all the 
A n are zero and the b is allowed to vary with time, the transient solution will pass at 
each instant through a steady-state configuration corresponding to the instantaneous b. 
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This is a quasi-steady solution. What is required is to find the variation of b with time 
for a given set of imposed conditions. Then figures 8 and 9 can be used to obtain the 
instantaneous quasi-steady configuration and heat conducted through the frozen region. 

The time variation of b is found from equation (57) by letting all the A n = 0. This 
gives 



OO 


-tf(l + q) ^ ^n^n+p 

n=0 



Since only q!q and (Sq are nonzero for A R = 0, their values from equations (62) and (63) 
are substituted to give 


b(-Ab) E T (l) _ 2?rAb 2 n 

7= o °> °>° T= / / A 

InVl - b 2 1 " b In Vl - b 2 K;(Yl-b 2 j 


( 66 ) 


From equation (56) 


^0. «- EG 1 )^ 1 [<-«»- a 

j=o 


or 


j(l) _ h 2 H^ + (2 - b 2 l 

J 0 , 0,0 " b - 1,0 + V D ' H 0 , 0 + W 0,0 


The H T s are found from equation (55) giving 


0 - b22K -l + < 2 - b2 > 2K 0 + p 2K ( 

i„i / Tb 2 


Substituting Q and E from equation (38) into equation (66) gives 
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By virtue of equation (53) the K n are only functions of b. Also K_j is seen to equal 
K^. Then the variables in equation (67) can be separated to give 




is the complete elliptic integral of the second kind. Also 
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du> 


i 


1 - b 2 sin 2 oo 


= K(b) 


These relations are substituted into equation (68). The resulting equation is simplified 
and then integrated from 0 = 0 to 0 to give 



K (l/l - k 2 j 

E(k)lnVl - k 2 + k 2 K(k) 

dk - 17 

(in/l-k 2 ) 2 

kA - k 2 j+ In Vl - k 2 

2A 2 


0 


(69) 


The dummy variable k has been introduced on the left to avoid confusion with the run- 
ning variable b. 

The integration starts at corresponding to any initial profile as shown in 

figure 8 resulting from imposed conditions A”. At time zero the imposed conditions are 
changed from A” to A. The denominator of the integrand goes to zero and hence 
0 — oo when b is such that 


A = 


In 


Vl - b 2 


bK^Vl - b 2 


^ initial * 


But this is exactly the condition for the steady-state layer corresponding to A as given 
by equation (64). When the cold plate is initially at t^ and there is no frozen layer, then 

lai ' rb 

If in equation (69) the integral is carried out from a lower limit of unity J , then 

then the integral for any initial b can be found by subtraction, that is, * 

_ y*^ _ J ' b initial 


/ 


^initial 


Hence by having integrated from 1 to arbitrary values of 


b, the quasi- steady time can be found for a layer to grow from any initial state to another 
state. The growth rate will of course be a function of A. 
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Results from Quasi-Steady Solution 


Figure 10 shows the variation of b as a function of dimensionless time for the 
quasi-steady solution. Curves are given for various values of A that are imposed at 
the beginning of the transient and maintained throughout the transient growth. As time 
proceeds, the b reaches a constant value corresponding to the steady-state frozen pro- 
file for the particular value of A that is imposed. 

The sets of curves in figures 7 to 10 can be used as follows to obtain the quasi- 
steady solution. First consider the case where initially the wall is at the freezing tem- 
perature t = t f . There will be no frozen layer on the plate and A - = °° since t f - t,„ 
is in the denominator of A. When A” — figure 7 shows that b — 1. This is consist- 
ent with figure 8 which shows that for b = 1 there is no frozen layer on the plate. To 
initiate the transient, at time 0 = 0 + the physical conditions are changed so that 


ha H " t f 


k ^ - ‘w 


goes from A = °° to the constant value of A that is maintained throughout the tran- 
sient. For that particular A the variation of b with 0 is found from figure 10. Then 
the time dependent profiles and heat flow during the transient are found from figures 8 
and 9. 

Now consider the case where A" is finite, that is, there is initially a steady-state 
layer on the plate. Corresponding to the value A” there is a b. n .j. . ^ that can be found 
from figure 7, and an initial frozen profile from figure 8. To start the transient, a new 
value of A is imposed and the transient variation of b will follow the curve in figure 10 
corresponding to the imposed A. To find the value of b after a given time lapse, the 
abscissa in figure 10 must be interpreted as being the actual 0 that has passed since 
the initiation of the transient plus the value of 0 corresponding to b^^.^ as obtained 
on the curve for A in figure 10. With the b known as a function of time from the onset 
of the transient, the heat flow and profile shape can be found as a function of time from 
figures 9 and 8. 
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Dimensionless time, 0 = .A- ~ 
kpA tj 

(b) Range of A, 0. 7 to 1.9. 
Figure 10. - Continued. 
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COMPUTATION OF TRANSIENT SOLIDIFICATION 

Having examined the steady and quasi-steady solutions, the general results will now 
be considered. The final result of the transient analysis is given by equation (57). This 
represents a set of equations, each equation corresponding to a different value of p. 

The mapping parameter b and the coefficients A n are unknown functions of time which 

must be found to evaluate the profiles given by equations (58). Since the set of equa- 

• • 

tions (57) involve the time derivatives b and (n = 0, 1, 2, . . . ), these are simulta- 
neous first-order ordinary differential equations. Equations (58) for the X_ and Y e 
coordinates of the interface each contain an infinite series. When the A are evaluated, 
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they decrease in size as n is increased so that the series can be terminated after a fi- 
nite number of terms. 

To find the b(O) and A n (0), the set of equations (57) is truncated at a value of p 
equal to 1 plus the number of A n desired. For example, if only b and Aq, . . . , Aq 
are retained, the only equations of the set which are retained are those for which 

QQ 

p = 0, 1, . . ., 4. The series ^ would extend only to n = 3. This set of five equa- 

n=0 

tions is solved using the Runge-Kutta method for simultaneous first-order differential 
equations. Since the series in A n has been truncated, the values of b and A n may 
not be accurate. A larger number of terms is then included, say Aq, . . . , Ag and 
p = 0, 1, . . ., 9. The first four A n are then checked against the A n found by using 
only five simultaneous equations. This procedure is continued until increasing the num- 
ber of equations does not change the values of b(0) and A R (0) . Also there must be a 
large enough number of A n so that the series in equations (58) have converged. 

The numerical computer program that was utilized could solve 20 simultaneous dif- 
ferential equations, so the present results are limited to cases that converged with 19 or 
less A . It was found that a larger number of A n would be required when the transient 
solidification started from a thin initial layer. It will be shown subsequently, however, 
that the quasi- steady solution works very well for thin initial layers so that the present 
theory does provide a means for computing transients starting from an initial condition 
of any steady-state frozen region or a bare plate. 

Now consider in more detail how equation (57) is written out for each p value. All 
of the quantities appearing in these equations such as a, E, and J must be replaced by 
their appropriate expressions in terms of b and A n . The expression for E is given 
by equation (38), and the expressions for the a n by equation (37). With the a n known, 
equation (49) is used to find each of the j3 n . Each /3 can be found successively by 
writing out equation (49) for successive values of n starting with n = 0. For example 

A n A a A A A A 

with n = 0, /3 q = cTq; then for n = 1, PqPi + PiPq = so that having found /3 q the 
can be found. Continuing, letting n = 2 gives an equation that can be solved for /3„ 

@ 0^2 + ^ 1 % + ^ 2^0 = a 2 

The „ given by equation (56) are only functions of b. They depend on the H 

n, k, p 

functions that are defined in equations (55). These depend on the K n which are the defi- 
nite integrals defined in equation (53). It was necessary to develop a table of K n values 
for various n and b by carrying out the definite integrals numerically. 

To obtain a transient solidification solution, initial conditions must be specified. The 
transients computed here are all started from an initial steady-state layer. The magni- 
tude of the physical quantity 
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that establishes this initial layer is called A - . From figure 7 the value of A” corre- 
sponds to a value of b which will be termed b- n ^ a1 . The initial shape of the solidified 
region corresponding to this value of b. n . f . a1 can be found from figure 8. At time 0 + 
a new physical condition specified by the value of A, is imposed and this remains fixed 
throughout the transient. For an initial steady-state layer, all the A r = 0 at 0 = 0 + . 
Hence the solution of equation (57) proceeds from the following initial conditions: 


b = t>i n itiaJL ( corres P on ding to initial steady profile) 

A n = ° 


ha h - *f _ / A ~ 0<O 

ktf-tw l A 0 > 0 


As time becomes large, the solidified layer approaches a steady state. The steady- 
state value of b is completely determined by the value of A from the curves of fig- 
ure 7, and the steady-state profile can be found at this value of b by use of figure 8. 

When the variations with time of b and the A n are known, equations (58) are used 
to compute the transient shapes of the ice layer. The transient heat flow through the 
layer is found from equation (59). 


RESULTS AND DISCUSSION OF TRANSIENT SOLUTION 
Solidification Computed from Transient Solution 

Several groups of transient growth curves are shown in figures 11 and 12 to demon- 
strate the nature of the transient solution. In figure 11, before the transients begin 
there is an initial steady-state layer that has been established by having the cooling pa- 
rameter equal to A” = 1.47 (corresponding to the initial condition = 0. 995). 

Then the cooling parameter is suddenly changed to another value, this being A = 0. 5 for 
figure 11(a). The change in A could physically correspond to changing the liquid heat 
transfer coefficient, the liquid temperature, or the temperature of the cooled plate. The 
frozen layer then grows and the profiles are shown for various values of the dimension- 
less time until a new steady state is achieved. On figures 11(b) and (c) the values of A 
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Dimensionless coordinate, y/a 






Dimensionless coordinate, y/a 





Dimensionless coordinate, y/a 




0 .5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

Dimensionless coordinate, x/a 


la) b initial = 0-98; A" ■ 1.04. 

Figure 12. - Transient solidification starting from three different initial layers and going to same final layer; A - 0. 1 
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(c) b| njt j a | ■ 0.8; A' ■ 0.366. 
Figure 12. - Concluded. 


are respectively 0. 3 and 0.2. A smaller value of A corresponds to increased cooling 
(increased t^ - t w ) or a decreased convection of energy h(t^ - t^) to the frozen interface. 
Consequently, for smaller A the steady- state layer (at © — °°) is larger. 

In figure 12 the final cooling parameter A is equal to 0. 1 for all three cases so 
that the steady-state layers are all the same. The transients begin, however, from 
various initial layers corresponding to cooling parameters prior to the transient A" of 
1.04, 0.560, and 0.366. 

The results in figures 11 and 12 give the reader a quantitative idea of the rate at 
which the solidification occurs. The rate is of course most rapid at early times when 
the frozen region has the least thickness and hence the least resistance to heat flowing 
through it. As the frozen region becomes large compared with the width of the cooled 
plate (x/a = 1) the frozen shape becomes circular, tending toward the axisymmetric solu- 
tion where the heat removal would be at a line sink at the center of the solidified region. 
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Figure 13 shows the heat flow through the solidified layer and into the wall for the 
six transients in figures 11 and 12. In figure 13(a) all the initial layers are the same 
and, hence, Q/2k(tj - t w ) starts from the same value. A small value of A corresponds 
to a large final solidified region and, consequently, a low steady- state heat flow, hi fig- 
ure 13(b), corresponding to the transients in figure 12, the A is the same for all cases 
and, hence, the same steady-state heat flow is reached. A large A” corresponds to a 
thin initial layer and consequently a large heat flow at the beginning of the transient. 



Figure 13. - Transient heat flow through frozen layer. 
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Comparison With Quasi-Steady Solution 


Figure 14 shows three different transient solutions and compares the layer shapes 
with those computed from the quasi-steady solution. These cases had the largest devia- 
tions that were obtained between the transient and quasi-steady solutions. It is evident 
that the transient profiles are always quite close to those predicted by the quasi-steady 
analysis. 

As mentioned previously, the series in the transient solution converged less rapidly 
for thin initial layers, and hence transient solutions were not carried out for thinner ini- 
tial layers than those shown in figure 11. However, thin layers as shown by figure 8 
would be almost one-dimensional except for the curved portion of the interface near the 
edge of the plate x/a ~ 1. Hence, when starting a transient from a thin layer, the tran- 
sient profiles must pass through a series of instantaneous steady states as both the tran- 
sient and quasi- steady layers are very flat. Consequently, the quasi- steady solution 



(a) Thin initial layer; bj n j tia j = 0.995; small cooling; A = 0.5. 
Figure 14. - Comparison of transient and quasi-steady solutions. 
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(b) Thin initial 


— Transient solution 

Quasi-steady solution 



■; bf nit j a | = 0.995; moderate cooling; A = 0.2. 
Figure 14. - Continued. 






(c) Thin initial layer; bj n | tia | = 0. 98; large cooling; A = 0. 1. 

Figure 14. - Concluded. 

should be a good engineering approximation for all cases. Thus for engineering calcula- 
tions it is a simple procedure to obtain transient results by use of figures 7 to 10 as de- 
scribed in the section entitled QUASI-STEADY SOLUTION. 


CONCLUSIONS 

A conformal mapping method was developed and applied to a particular case of tran- 
sient two-dimensional solidification. The configuration in the physical plane is obtained 
by a quadrature that involves the properties of the solidified region in a potential plane 
and in a temperature derivative plane. To perform the integration, the configurations 
in these planes are both mapped conformally into a fixed region of an intermediate plane. 
Since the regions in the potential and derivative planes depend on time, the mapping 
functions also contain time varying quantities. These quantities are found by satisfying 
the physical conditions at the moving interface. 
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A quasi-steady solution was also carried out. In this solution the profiles pass 
through a series of instantaneous steady-state shapes. The transient results from the 
general analysis were found to be predicted within engineering accuracy by the quasi- 
steady solution. The quasi-steady solution is provided in graphical form and can be used 
to quickly predict the solidification behavior. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, August 28, 1969, 

129-01. 
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APPENDIX - COMPUTATION OF THE INTEGRAL FOR 1 ( 0 ., 0) 


To obtain the expression for Z (O, 0) given by equation (25), it is necessary to eval- 


uate 


Z - A = 


°Q _ 

I'/ 




n , 

r 


l 


n=0 


(i +y ) z -(i-b^)(i- r r 


(Al) 


To this end, examine the form of the integral 


/ 


n , 

r 

Vi 


(A2) 


where 


Vx = V(1 +y) 2 - 


(l + y) “ (1 - b )(1 - y) 


for a few values of n. Any standard integral table shows that the results of these inte- 
grations are as follows: 


— ln ^2 - b 2 + b 2 y + b Vx ) + const for n = 0 


(A3) 


Vx _ 4_ - 2 b^ ln - b 2 + b 2 y + bVx) + const 


for n = 1 


(A4) 


b 2b‘ 


2b 2 y - 6(2 - b 2 ) 

4b 4 


8b 


Vx + 3(4 - 2b 2 ) 2 ^4 b^ ln - b 2 + b 2 y + bVx) + const for n = 2 


(A5) 


These integrations indicate that the n integral can be expressed in the form 
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/ 


n 


n-1 


3— dy = Vx ^ + C n In (2 - b 2 + b 2 y + bVx) + const r > 0 

Vx r=0 


(A 6) 


where the /3^ and C n are coefficients that will now be found. Differentiate equa- 
tion (A6) with respect to y (note that X contains y) 


4 - V5 £ + - « - 2» 2 g £ 

VX r=l 2VX r=0 


n 


2 - b 2 + b 2 y + bVx 

Multiply through by Vx and simplify the last term to obtain 


b 2 + - JL- (2b 2 y + 4 - 2b 2 ) 


\Vx 


y n = X £ r^y- 1 + (b 2 y + 2 - b 2 ) £ £y r + C n b 
r=l r=0 

Substitute for X = (1 + y) 2 - (1 - b 2 )(l - y) 2 = b 2 y 2 + (4 - 2b 2 )y + b 2 to give 

y 11 = [b 2 y 2 + (4 - 2b 2 )y + b 2 ] £ r j^y r_1 + (b 2 y + 2 - b 2 ) £ $y r + C n b 

r=l r=0 

= b 2 2 (1 + r)/£y r+1 + (2 - b 2 ) £ (1 + 2r)^y r + b 2 £ + C n b 

r=0 r=0 r=0 


(A7) 


Terms in the same power of y will be equated to determine the jEjJJ and the C . To 
this end equation (A7) is rearranged further to obtain 


n 


y” = b 2 £ r^_ lV r + (2 - b 2 )^g + (2 - b 2 ) £ (1 + 2r)^y r + b 2 £ r^y 1 " 1 + C n b 
r=l r=l r=l 
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which can be put in the form 


y 11 = (2 - b 2 )/?Q + b 2 n^_ ir n + b 2 (n - l)^_ 2 y n_1 + (2 - b 2 )[l + 2(n - DJ^y 


n-1 


n-2 


n-2 


+ 


S [ b2 <-i * < 2 - h 2 )* 1 + + b 2 E (r + l)Cir r * c n b 


r=l 


r=0 


y 11 = b 2 n^_ir n + [b 2 (n - 1)£_ 2 + (2 - b 2 )(2n - lJjSj.Jy 11 " 1 


n-2 


+ ^ [b 2 r^_j + (2 - b 2 )(l + 2r)/£ + b 2 (r + l)/3j +1 ]y r + (2 - b 2 )(§ + b 2 ^ + C n b 
r=l 

(A 8) 


Equating the terms in y® gives the C n for n > 0 


(2 - b 2 )$ + b 2 # 


C n = - 


n = 1, 2, 3, . . . 


Comparing equations (A6) and (A3) gives the Cq as 


C - — 

0 b 


Equating the coefficients of y 11 in equation (A8) gives 


£ n , - 1 
'n-1 


b 2 n 


x 

From the (n - l) in order terms 


b 2 (n - l)/3g_ 2 + (2 - b 2 )(2n - l)^_j = 0 


which gives 


68 



I 


0 


n 

n-2 


2 - b 2 (2n - 1) 
b 4 n(n - 1) 


For 1 < r < n - 2 

b 2 ^^ + (2 - b 2 )(l + 2r)/3” + b 2 (r + 1)0^ +1 = 0 

which gives 


0 


n 

r-1 " 




r + 1 
r 


P 


n 

r+1 


This shows the way in which the coefficients in equations (26) and (27) were obtained. 
To evaluate equation (Al), insert equation (A 6) to obtain 


Z - A = 



a. 


n 


n-l 

Vx Yj 0r yP + C n ln ( 2 ' b2 + b2 r + b Vx) 
r=0 


-.ST 


n=0 


y =0 


= In 


~2 - b 2 + b 2 ^ 2 + bVx(^) 


OO oo n- 1 OO 

E “„c n+ E % E *? n2r - b E v 

n=0 n=l r=0 n=l 


(A9) 


This is the form of equation (25a). 
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